Participant Fluctuation Amplitude Data

fluctuations.glasser.pnc <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/PNC/glasserfluctuations_demographics_finalsample.csv") #fMRI + demographics, generated with fitGAMs_fluctuationamplitude_age.R

Glasser Labels

glasser.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/parcellations/surface/glasser360_regionlist.csv", header = T) #glasser region list

Sensorimotor-Association Axis

S.A.axis.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo, vertex-wise axis average ranks
S.A.axis <- as.data.frame(cbind(rank(S.A.axis.cifti$data), names(S.A.axis.cifti$Parcel)))
colnames(S.A.axis) <- c("SA.axis","orig_parcelname")
S.A.axis <- merge(S.A.axis, glasser.parcel.labels, by="orig_parcelname", sort = F)
S.A.axis$SA.axis <- as.numeric(S.A.axis$SA.axis)
rm(S.A.axis.cifti)

Intracortical Myelination Developmental Data

myelin.partialR2.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Myelin/hcpd_n628_myelin_sAge_partial_bayes_r2.pscalar.nii") #T1/T2 ratio - age effect size
myelin.maxdev.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Myelin/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin <- as.data.frame(cbind(myelin.partialR2.cifti$data, myelin.maxdev.cifti$data))
colnames(myelin) <- c("myelin.partialR2","myelin.maxdev")
myelin$label <- glasser.parcel.labels$label
rm(myelin.partialR2.cifti)
rm(myelin.maxdev.cifti)

SNR Mask

SNR.mask <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/parcellations/surface/SNRmask_glasser360.csv")

Spin Test Parcel Rotation Matrix

source("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/R/rotate.parcellation.R")
source("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/R/perm.sphere.p.R")
glasser.coords <- read.table("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/sphere_HCP.txt", header=F) #coordinates of glasser parcel centroids on the freesurfer sphere
perm.id.full <- rotate.parcellation(coord.l = as.matrix(glasser.coords[1:180,]), coord.r = as.matrix(glasser.coords[181:360,]), nrot = 10000) #rotate the glasser parcellation 10,000 times on the freesurfer sphere to generate spatial nulls for spin-based permutation significance testing 
saveRDS(perm.id.full, "/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds")
perm.id.full <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins

GAM Results

gam.age.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_statistics_glasser.csv") #GAM age smooth statistics, generated with fitGAMs_fluctuationamplitude_age.R

#create df.dev with GAM statistics, S-A axis, plasticity measures, and SNR mask
df.list <- list(gam.age.glasser, S.A.axis, myelin, SNR.mask) #dfs to merge
df.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
df.dev <- df.dev %>% select(label, orig_parcelname, everything()) #order columns
cols = c(3:15)    
df.dev[,cols] = apply(df.dev[,cols], 2, function(x) as.numeric(as.character(x))) #format numerics

#create df.dev.spin formatted for spatial permutation testing
df.dev.spin <- rbind(df.dev[181:360,], df.dev[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.dev.spin$GAM.age.partialR2[df.dev.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.dev.spin$minage.decrease[df.dev.spin$SNR.mask == 0] <- NA 

df.dev <- df.dev %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analyses
df.dev$SA.axis.bin <- as.numeric(cut2(df.dev$SA.axis, g=10)) #divide the S-A axis into 10 ranked bins
gam.fitted.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_predictedfits_glasser.csv") #GAM predicted fits, generated with fitGAMs_fluctuationamplitude_age.R
gam.fitted.glasser <- inner_join(gam.fitted.glasser, df.dev, by="label", sort = F)
gam.smoothestimates.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_smoothestimates_glasser.csv") #GAM smooth estimates, generated with fitGAMs_fluctuationamplitude_age.R
gam.smoothestimates.glasser <- inner_join(gam.smoothestimates.glasser, df.dev, by="label", sort = F)
gam.derivatives.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_derivatives_glasser.csv") #GAM true model age smooth derivatives, generated with fitGAMs_fluctuationamplitude_age.R
gam.derivatives.glasser <- left_join(gam.derivatives.glasser, S.A.axis, by="label", sort = F)
gam.derivatives.glasser <- left_join(gam.derivatives.glasser, SNR.mask, by="label", sort = F)
gam.derivatives.glasser <- gam.derivatives.glasser %>% filter(SNR.mask != 0) 
gam.sex.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_agebysex_interaction_glasser.csv") #GAM age spline by sex interaction statistics, generated with fitGAMs_fluctuationamplitude_environment.R

#create df.sex with GAM results, S-A axis, and SNR mask
df.list <- list(gam.sex.glasser, S.A.axis, SNR.mask) #dfs to merge
df.sex <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
gam.env.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_environment_statistics_glasser.csv") #GAM environment statistics, generated with fitGAMs_fluctuationamplitude_environment.R

#create df.env with GAM results, S-A axis, and SNR mask
df.list <- list(gam.env.glasser, S.A.axis, SNR.mask) #dfs to merge
df.env <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 

#create df.env.spin formatted for spatial permutation testing
df.env.spin <- rbind(df.env[181:360,], df.env[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.env.spin$GAM.env.tvalue[df.env.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.env.spin$GAM.env.edu.tvalue[df.env.spin$SNR.mask == 0] <- NA

df.env <- df.env %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analyses

Age-dependent changes in intrinsic fMRI fluctuations vary across the cortex

Age spline by sex interactions

Number of cortical regions show significance age*sex interactions

sum(p.adjust(df.sex$GAM.agexsex.pvalue, method=c("fdr")) < 0.05) #FDR-corrected p-values age by sex interaction testing whether the age spline differs between males and females
## [1] 0

Cortical smooth functions

gam.smoothestimates.glasser.lh <- gam.smoothestimates.glasser[33401:67200,]
ggplot(gam.smoothestimates.glasser.lh,aes(age,est,group=index,color=GAM.age.partialR2)) + 
  geom_line(size=.3, alpha = .8) + 
  paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme_classic() +
  labs(x = "\nAge", y = "Fluctuation Amplitude (zero centered)\n" ) +
  theme(legend.position = "none") +
  theme(axis.text = element_text(size=6, family = "Arial", color = c("black")), axis.title = element_text(size=6, family = "Arial", color = c("black")), axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure1/Regional_agesmooths.jpg", device = "pdf", dpi = 500, width = 2.55 , height = 2.3)

Region-specific developmental trajectories and derivatives

Region V1

V1.smooth <- gam.fitted.glasser %>% filter(label == "rh_R_V1")
V1.derivatives <- gam.derivatives.glasser %>% filter(label == "rh_R_V1")

scatterplot <- ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = rh_R_V1)) +
  geom_point(color="#FBDC9D", size = .3) +
  geom_line(data = V1.smooth, aes(x = age, y = fitted), color="#F8B57B", size = 1) +
  geom_ribbon(data = V1.smooth, aes(x = age, y = fitted,  ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#F8B87D",) +
  labs(x="\nAge", y="Fluctuation Amplitude\n") +
  theme_classic() +
  theme(
  axis.text = element_text(size=6, family = "Arial", color = c("black")),
  axis.title.x=element_text(size=6, family ="Arial", color = c("black")),
  axis.title.y=element_text(size=6, family ="Arial", color = c("black")),
  axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(breaks=c(1.0, 2.0), labels = c(sprintf("%.1f", round(1, 2)), sprintf("%.1f", round(2, 2))))

derivplot <- ggplot(data = V1.derivatives) + 
  geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
  scale_fill_gradient(low = "#EFAF77", high = "#FFE0A1", limits = c(-0.0323371,-.01), na.value = "white") +
  labs(x="\nAge", fill = "Derivative") + 
  theme_classic() +
  theme(legend.position = "none") +
  theme(axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size=6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

allplots <- list(scatterplot,derivplot)
V1.plot <- cowplot::plot_grid(rel_heights = c(16,3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
V1.plot 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure1/V1_smooth.jpg", device = "pdf", dpi = 500, width = 2.1 , height = 1.85)

Region p24pr

p24pr.smooth <- gam.fitted.glasser %>%filter(label == "rh_R_p24pr")
p24pr.derivatives <- gam.derivatives.glasser %>% filter(label == "rh_R_p24pr")

scatterplot <- ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = rh_R_p24pr)) +
  geom_point(color = "#e06d65", size = .3) +
  geom_line(data = p24pr.smooth, aes(x = age, y = fitted), color="#D54D55", size = 1) +
  geom_ribbon(data = p24pr.smooth, aes(x = age, y = fitted,  ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#D54D55",) +
  labs(x="\nAge", y="Fluctuation Amplitude\n") +
  theme_classic() +
  theme(
  axis.text = element_text(size=6, family = "Arial", color = c("black")),
  axis.title.x=element_text(size=6, family ="Arial", color = c("black")),
  axis.title.y=element_text(size=6, family ="Arial", color = c("black")),
  axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(breaks=c(0.5, 1.0))

derivplot <- ggplot(data = p24pr.derivatives) + 
  geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
  scale_fill_gradient(high = alpha("#DF5E54",0.7), low = "#D54D55", na.value = "white", limits = c(min(p24pr.derivatives$significant.derivative),-.001)) +
  labs(x="\nAge", fill = "Derivative") + 
  theme_classic() +
  theme(legend.position = "none") +
  theme(axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size=6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

allplots <- list(scatterplot,derivplot)
p24pr.plot <- cowplot::plot_grid(rel_heights = c(16,3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
p24pr.plot 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure1/p24pr_smooth.jpg", device = "pdf", dpi = 500, width = 2.1 , height = 1.85)

Region IFSa

IFSa.smooth <- gam.fitted.glasser %>%filter(label == "lh_L_IFSa")
IFSa.derivatives <- gam.derivatives.glasser %>% filter(label == "lh_L_IFSa")

scatterplot <- ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = lh_L_IFSa)) +
  geom_point(color="#893774", size = .3) +
  geom_line(data = IFSa.smooth, aes(x = age, y = fitted), color="#381043", size = 1) +
  geom_ribbon(data = IFSa.smooth, aes(x = age, y = fitted,  ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#381043",) +
  labs(x="\nAge", y="Fluctuation Amplitude\n") +
  theme_classic() +
  theme(
  axis.text = element_text(size=6, family = "Arial", color = c("black")),
  axis.title.x=element_text(size=6, family ="Arial", color = c("black")),
  axis.title.y=element_text(size=6, family ="Arial", color = c("black")),
  axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(breaks=c(1.0, 2.0), labels = c(sprintf("%.1f", round(1, 2)), sprintf("%.1f", round(2, 2))))

derivplot <- ggplot(data = IFSa.derivatives) + 
  geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
  scale_fill_gradient2(high = "#6E195E", low = "#cf93c4", midpoint = 0, mid = "white", na.value = "white") +
  labs(x="\nAge", fill = "Derivative") + 
  theme_classic() +
  theme(legend.position = "none") +
  theme(axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size=6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

allplots <- list(scatterplot,derivplot)
IFSa.plot <- cowplot::plot_grid(rel_heights = c(16,3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
IFSa.plot 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure1/IFSa_smooth.jpg", device = "pdf", dpi = 500, width = 2.1 , height = 1.85)

The development of intrinsic fMRI fluctuations is linked to the maturation of a key regulator of plasticity

Spatial development

Fluctuation amplitude age effects map

ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2), position = c("stacked"), hemisphere = c("right")) + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure2/Brainmap_function.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

T1/T2 ratio age effects map

ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=myelin.partialR2), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = 1, limits = c(0.0,0.3), oob = squish) +
  theme(legend.text = element_blank())
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##    atlas type  hemi  side  region label    orig_parcelname GAM.age.Fvalue
##    <chr> <chr> <chr> <chr> <chr>  <chr>    <chr>                    <dbl>
##  1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V1  L_V1_ROI                 27.6 
##  2 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_MST L_MST_ROI                14.7 
##  3 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V6  L_V6_ROI                  7.44
##  4 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V2  L_V2_ROI                 33.4 
##  5 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V3  L_V3_ROI                 27.8 
##  6 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V4  L_V4_ROI                 45.6 
##  7 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V8  L_V8_ROI                 29.3 
##  8 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_4   L_4_ROI                  20.6 
##  9 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_3b  L_3b_ROI                 16.2 
## 10 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_FEF L_FEF_ROI                13.3 
## # … with 159 more rows, and 13 more variables: GAM.age.pvalue <dbl>,
## #   GAM.age.partialR2 <dbl>, Anova.age.pvalue <dbl>, age.onsetchange <dbl>,
## #   age.peakchange <dbl>, minage.decrease <dbl>, maxage.increase <dbl>,
## #   age.maturation <dbl>, SA.axis <dbl>, myelin.partialR2 <dbl>,
## #   myelin.maxdev <dbl>, SNR.mask <dbl>, SA.axis.bin <dbl>

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure2/Brainmap_myelin.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

Correlation (rho + p.spin) between fluctuation amplitude age effects and T1/T2 ratio age effects

cor.test(df.dev$GAM.age.partialR2, df.dev$myelin.partialR2, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$GAM.age.partialR2 and df.dev$myelin.partialR2
## S = 10557746, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.6699693
perm.sphere.p(df.dev.spin$GAM.age.partialR2, df.dev.spin$myelin.partialR2, perm.id.full,corr.type='spearman')
## [1] 0.00045

Correlation plot

ggplot(df.dev, aes(x = myelin.partialR2, y = GAM.age.partialR2, fill = GAM.age.partialR2)) + 
geom_point(aes(color = GAM.age.partialR2), shape = 21, size = 1.5) +
paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
labs(x="\nMyelin Age Effect", y="Fluctuation Amplitude Age Effect\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .75) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22))
## `geom_smooth()` using formula 'y ~ x'

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure2/Ageeffects_correlationplot.jpg", device = "pdf", dpi = 500, width = 2.68 , height = 2.175)
## `geom_smooth()` using formula 'y ~ x'

Temporal development

Fluctuation amplitude age of decrease onset

ggseg(.data = df.dev, atlas = "glasser", mapping = aes(fill = minage.decrease), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill = "") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10, 17.5), oob = squish) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure2/Brainmap_function_agedecrease.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

T1/T2 ratio age of maximal increase

ggseg(.data = df.dev, atlas = "glasser", mapping = aes(fill = myelin.maxdev), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill = "") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10,17.5), oob = squish) +
  theme(legend.text = element_blank())
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##    atlas type  hemi  side  region label    orig_parcelname GAM.age.Fvalue
##    <chr> <chr> <chr> <chr> <chr>  <chr>    <chr>                    <dbl>
##  1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V1  L_V1_ROI                 27.6 
##  2 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_MST L_MST_ROI                14.7 
##  3 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V6  L_V6_ROI                  7.44
##  4 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V2  L_V2_ROI                 33.4 
##  5 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V3  L_V3_ROI                 27.8 
##  6 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V4  L_V4_ROI                 45.6 
##  7 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V8  L_V8_ROI                 29.3 
##  8 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_4   L_4_ROI                  20.6 
##  9 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_3b  L_3b_ROI                 16.2 
## 10 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_FEF L_FEF_ROI                13.3 
## # … with 159 more rows, and 13 more variables: GAM.age.pvalue <dbl>,
## #   GAM.age.partialR2 <dbl>, Anova.age.pvalue <dbl>, age.onsetchange <dbl>,
## #   age.peakchange <dbl>, minage.decrease <dbl>, maxage.increase <dbl>,
## #   age.maturation <dbl>, SA.axis <dbl>, myelin.partialR2 <dbl>,
## #   myelin.maxdev <dbl>, SNR.mask <dbl>, SA.axis.bin <dbl>

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure2/Brainmap_myelin_ageincrease.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

Correlation (rho + p.spin) between fluctuation amplitude age of decrease onset and T1/T2 ratio max age of increase

#All regions
cor.test(df.dev$minage.decrease, df.dev$myelin.maxdev, method=c("spearman"))
## Warning in cor.test.default(df.dev$minage.decrease, df.dev$myelin.maxdev, :
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$minage.decrease and df.dev$myelin.maxdev
## S = 1677077, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6382729
perm.sphere.p(df.dev.spin$minage.decrease, df.dev.spin$myelin.maxdev, perm.id.full,corr.type='spearman')
## [1] 0.00125
#Age of decrease onset > 8.5
late.decrease <- df.dev %>% filter(minage.decrease > 8.5)
cor.test(late.decrease$minage.decrease, late.decrease$myelin.maxdev, method=c("spearman"))
## Warning in cor.test.default(late.decrease$minage.decrease,
## late.decrease$myelin.maxdev, : Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  late.decrease$minage.decrease and late.decrease$myelin.maxdev
## S = 349670, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6402458
late.decrease.spin <- df.dev.spin
late.decrease.spin$minage.decrease[late.decrease.spin$minage.decrease < 8.5] <- NA
perm.sphere.p(late.decrease.spin$minage.decrease, late.decrease.spin$myelin.maxdev, perm.id.full,corr.type='spearman')
## [1] 0.01565

Correlation plot

ggplot(late.decrease, aes(x = myelin.maxdev, y = minage.decrease, fill = minage.decrease)) + 
geom_point(aes(color = minage.decrease), shape = 21, size = 1.2) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10, 17.5), oob = squish) +
paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10, 17.5), oob = squish) +
labs(x="\nMyelin Age Effect", y="Fluctuation Amplitude Age Effect\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .75) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22)) +
scale_y_continuous(breaks=c(10, 12, 14, 16, 18, 20))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure2/Temporaleffects_correlationplot.jpg", device = "pdf", dpi = 500, width = 2.56 , height = 2.175)

The development of intrinsic fMRI fluctuations aligns with the hierarchical sensorimotor-association axis

Sensorimotor-association axis

ggseg(.data = df.dev, atlas = "glasser", mapping = aes(fill = SA.axis), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill = "") +
  scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure3/SAaxis.jpg", device = "pdf", dpi = 500, width = 3.65, height = 2.95)

Age effects along the S-A Axis

Correlation (rho + p.spin) between fluctuation amplitude age effects and the sensorimotor-association axis

cor.test(df.dev$GAM.age.partialR2, df.dev$SA.axis, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$GAM.age.partialR2 and df.dev$SA.axis
## S = 2934378, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5358554
perm.sphere.p(df.dev.spin$GAM.age.partialR2, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')
## [1] 0.00215

Temporal change along the S-A Axis

Correlation (rho + p.spin) between fluctuation amplitude age of decrease onset and the sensorimotor-association axis

cor.test(df.dev$minage.decrease, df.dev$SA.axis, method=c("spearman"))
## Warning in cor.test.default(df.dev$minage.decrease, df.dev$SA.axis, method =
## c("spearman")): Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$minage.decrease and df.dev$SA.axis
## S = 2477911, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4655417
perm.sphere.p(df.dev.spin$minage.decrease, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')
## [1] 0.0388
cor.test(late.decrease$minage.decrease, late.decrease$SA.axis, method=c("spearman"))
## Warning in cor.test.default(late.decrease$minage.decrease,
## late.decrease$SA.axis, : Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  late.decrease$minage.decrease and late.decrease$SA.axis
## S = 314027, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6769174
perm.sphere.p(late.decrease.spin$minage.decrease, late.decrease.spin$SA.axis, perm.id.full,corr.type='spearman')
## [1] 0.0011

Sensorimotor-association axis binned GAM smooths

gam.smoothestimates.glasser$age <- round(gam.smoothestimates.glasser$age, 3)

meansmooth.bins <- gam.smoothestimates.glasser %>% group_by(age, SA.axis.bin) %>% do(est.mean = mean(.$est)) %>% unnest(cols = c(est.mean)) #calculate the average smooth estimate at each age within each of the 10 S-A axis bins
ggplot(data = meansmooth.bins, aes(x = age, y = est.mean, group = SA.axis.bin, color = as.factor(SA.axis.bin))) +
  geom_line(size = 1.5) +
  scale_color_manual(values = c("#FFC228", "#FFCA4E", "#FFD16A", "#FFDA89", "#FFE6B0", "#D7BCDA", "#BE93C4", "#9859A4", "#863E95", "#6F1382")) +
  theme_classic() +
  labs(x = "\nAge", y = "Fluctuation Amplitude (zero centered)\n", color = "Sensorimotor-Association \nAxis Percentile Bin") +
  theme(legend.position = c(.585, 0.94), legend.direction = "horizontal", legend.text = element_text(size = 6), legend.key.size = unit(.3, "cm"), legend.title = element_text(size = 6, family = "Arial", color = c("black"))) +
  theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size=.22), axis.ticks = element_line(size = .22)) + scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure3/SAaxis_bin_smooths.jpg", device = "pdf", dpi = 500, width = 3.34 , height = 2.15)

Principal developmental gradient

Principal component analysis to identify the spatial axis that captures the most variance in regional fluctuation amplitude developmental trajectories

gam.smoothestimates.glasser <- gam.smoothestimates.glasser %>% select(label, age, est)
gam.smoothestimates.long <- gam.smoothestimates.glasser %>% pivot_wider(names_from="label",values_from="est")
gam.smoothestimates.long <- gam.smoothestimates.long %>% select(-age) #ages by regions matrix for PCA

smooths.pca <- prcomp(gam.smoothestimates.long) #PCA

Variance explained by PC1

summary(smooths.pca)$importance[2,1]
## [1] 0.87026

Developmental PC1 correlation with the sensorimotor-association axis

PC1 <- as.data.frame(smooths.pca$rotation[,1]) #first principal component loadings
PC1$label <- row.names(PC1) 
colnames(PC1) <- c("PC1","label")
PC1$ranked <- rank(PC1$PC1)

cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  PC1$PC1 and df.dev$SA.axis
## S = 1896048, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.700093
df.dev.spin <- left_join(df.dev.spin, PC1, by="label")
perm.sphere.p(df.dev.spin$PC1, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')
## [1] 0

Correlation plot

PC1 <- left_join(PC1, S.A.axis, by="label")
ggplot(PC1, aes(x = SA.axis, y = PC1, fill = SA.axis)) + 
geom_point(aes(color = SA.axis), shape = 21, size = 1.2) +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis Rank", y="Principal Developmental Axis Loading\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .25) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure3/AxisPC1_correlationplot.jpg", device = "pdf", dpi = 500, width = 3.34 , height = 2.15)

Principal developmental axis brain map

ggseg(.data = PC1, atlas = "glasser", mapping=aes(fill = PC1), position = c("stacked"), hemisphere = "right") + 
theme_void() + 
labs(fill = "") +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.05, limits = c(-0.11,0.008), oob = squish) + theme(legend.text = element_blank())
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##    atlas type  hemi  side  region label        PC1 ranked orig_parcelname SA.axis
##    <chr> <chr> <chr> <chr> <chr>  <chr>      <dbl>  <dbl> <chr>             <dbl>
##  1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V1  -0.103      12 L_V1_ROI              1
##  2 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_MST -0.0635     92 L_MST_ROI            61
##  3 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V6  -0.0432    164 L_V6_ROI             13
##  4 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V2  -0.0940     18 L_V2_ROI              6
##  5 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V3  -0.0910     22 L_V3_ROI             18
##  6 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V4  -0.109       8 L_V4_ROI             33
##  7 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_V8  -0.0816     36 L_V8_ROI             34
##  8 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_4   -0.0699     72 L_4_ROI              16
##  9 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_3b  -0.0647     85 L_3b_ROI             10
## 10 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L_FEF -0.0383    181 L_FEF_ROI           141
## # … with 159 more rows

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure3/PC1.jpg", device = "pdf", dpi = 500, width = 3.65, height = 2.95)

Developmental PC1 correlation with anatomical, functional, and evolutionary cortical hierarchies

maps <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/S-A_ArchetypalAxis/Glasser360_MMP_WholeBrain/brainmaps_glasser.csv")
maps$label <- glasser.parcel.labels$label
maps <- merge(PC1, maps, by = "label")

Anatomical hierarchy

cor.test(maps$PC1, maps$T1T2ratio, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  maps$PC1 and maps$T1T2ratio
## S = 10161350, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.6072694
cocor.dep.groups.overlap(
  r.jk = (cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))$estimate), 
  r.jh = (abs(cor.test(maps$PC1, maps$T1T2ratio, method=c("spearman"))$estimate)), 
  r.kh = (abs(cor.test(maps$T1T2ratio, maps$SA.axis, method=c("spearman"))$estimate)), 
  n = 336, alternative = "two.sided", test = "hittner2003")
## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk = 0.7001 and r.jh = 0.6073
## Difference: r.jk - r.jh = 0.0928
## Related correlation: r.kh = 0.7785
## Group size: n = 336
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = 3.5208, p-value = 0.0004
##   Null hypothesis rejected

Functional hierarchy

cor.test(maps$PC1, maps$G1.fMRI, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  maps$PC1 and maps$G1.fMRI
## S = 2560212, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.595039
cocor.dep.groups.overlap(
  r.jk = (cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))$estimate), 
  r.jh = (abs(cor.test(maps$PC1, maps$G1.fMRI, method=c("spearman"))$estimate)), 
  r.kh = (abs(cor.test(maps$G1.fMRI, maps$SA.axis, method=c("spearman"))$estimate)), 
  n = 336, alternative = "two.sided", test = "hittner2003")
## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk = 0.7001 and r.jh = 0.595
## Difference: r.jk - r.jh = 0.1051
## Related correlation: r.kh = 0.8664
## Group size: n = 336
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = 5.0048, p-value = 0.0000
##   Null hypothesis rejected

Evolutionary hierarchy

cor.test(maps$PC1, maps$Evolution.Expansion, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  maps$PC1 and maps$Evolution.Expansion
## S = 4298202, p-value = 2.362e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3201328
cocor.dep.groups.overlap(
  r.jk = (cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))$estimate), 
  r.jh = (abs(cor.test(maps$PC1, maps$Evolution.Expansion, method=c("spearman"))$estimate)), 
  r.kh = (abs(cor.test(maps$Evolution.Expansion, maps$SA.axis, method=c("spearman"))$estimate)), 
  n = 336, alternative = "two.sided", test = "hittner2003")
## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk = 0.7001 and r.jh = 0.3201
## Difference: r.jk - r.jh = 0.38
## Related correlation: r.kh = 0.6026
## Group size: n = 336
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = 9.6478, p-value = 0.0000
##   Null hypothesis rejected

Neurodevelopment is hierarchical through adolescence

Age-specific fluctuation amplitude derivatives

Regional derivative by age plot

gam.derivatives.glasser.lh <- gam.derivatives.glasser[33401:67200,]
ggplot(gam.derivatives.glasser.lh,aes(age,derivative,group=index,color=SA.axis)) + 
  geom_line(size=1, alpha = .7) + 
  scale_color_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = 180) +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "\nAge", y = "Derivative\n", color = "Sensorimotor-Association Axis Bin") +
  theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22)) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

gam.derivatives.glasser$age <- round(gam.derivatives.glasser$age, 5)

Age 10 derivative map

age10.derivs <- gam.derivatives.glasser %>% filter(age == 10.03015)

ggseg(.data = age10.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Age10_derivatives.jpg", device = "pdf", dpi = 500, width = 2.95 , height = 1.95)

Age 10 derivative-axis correlation plot

ggplot(age10.derivs, aes(x = SA.axis, y = derivative, fill = derivative)) + 
geom_point(aes(color = derivative), shape = 21, size = .1) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
labs(x="\nSensorimtor-Association Axis Rank", y="Derivative\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .35) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_blank(), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .15), axis.ticks = element_line(size = .15)) +
scale_y_continuous(limits = c(-.05, .055), oob = squish)
## `geom_smooth()` using formula 'y ~ x'

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Age10_correlationplot.jpg", device = "pdf", dpi = 500, width = 1.45 , height = 1.25)
## `geom_smooth()` using formula 'y ~ x'

Age 15 derivative map

age15.derivs <- gam.derivatives.glasser %>% filter(age == 15.02429)

ggseg(.data = age15.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Age15_derivatives.jpg", device = "pdf", dpi = 500, width = 2.95 , height = 1.95)

Age 15 derivative-axis correlation plot

ggplot(age15.derivs, aes(x = SA.axis, y = derivative, fill = derivative)) + 
geom_point(aes(color = derivative), shape = 21, size = .1) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
labs(x="\nSensorimtor-Association Axis Rank", y="Derivative\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .35) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_blank(), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .15), axis.ticks = element_line(size = .15)) +
scale_y_continuous(limits = c(-.05, .055), oob = squish)
## `geom_smooth()` using formula 'y ~ x'

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Age15_correlationplot.jpg", device = "pdf", dpi = 500, width = 1.45 , height = 1.25)
## `geom_smooth()` using formula 'y ~ x'

Age 20 derivative map

age20.derivs <- gam.derivatives.glasser %>% filter(age == 20.01843)

ggseg(.data = age20.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025,0.005), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Age20_derivatives.jpg", device = "pdf", dpi = 500, width = 2.95 , height = 1.95)

Age 20 derivative-axis correlation plot

ggplot(age20.derivs, aes(x = SA.axis, y = derivative, fill = derivative)) + 
geom_point(aes(color = derivative), shape = 21, size = .1) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
labs(x="\nSensorimtor-Association Axis Rank", y="Derivative\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .35) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_blank(), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .15), axis.ticks = element_line(size = .15)) +
scale_y_continuous(limits = c(-.05, .055), oob = squish)
## `geom_smooth()` using formula 'y ~ x'

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Age20_correlationplot.jpg", device = "pdf", dpi = 500, width = 1.45 , height = 1.25)
## `geom_smooth()` using formula 'y ~ x'

Derivative by sensorimotor-association axis plot

derivative.colorbar <- c("#ff7e38", "#f29049", "#ffa25e", "#faab73", "#fabb73", "#ffcd94", "#ffd894", "#ffe6ab", "#ffe6ab", "#fff5de", "#ffffff", "#f79a97", "#f58682", "#EE6E69", "#C95665", "#BB4166", "#B72863", "#941348", "#7d0f65",  "#660c5f", "#4f0e4a")

ggplot(data=gam.derivatives.glasser,aes(x = age, y = SA.axis, group = index)) +
  geom_line(aes(color=derivative), size=1) +
  scale_color_gradientn(colours = derivative.colorbar, limits=c(-0.028, 0.028), oob = squish) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  labs(color="") +
   theme(axis.text = element_text(size = 7, family = "Arial", color = c("black")), axis.title = element_text(size = 7, family = "Arial", color = c("black")), axis.line = element_line(size =.22), axis.ticks = element_line(size = .22), legend.text = element_text(size = 6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))   

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Derivative_rank_age_plot.pdf", device = "pdf", dpi = 500, width = 3.7 , height = 4.9)

Age-specific correlation between derivates and the S-A Axis

Posterior derivative - sensorimotor-association axis correlations, N=10,000

deriv.SAaxis.posteriorcorrs <- read.table("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/SAaxis_posteriorderivative_correlation_byage_glasser.csv", header = F, sep = ",") #get the regional derivative-regional axis rank correlation at each age for all draws from the posterior distribution
colnames(deriv.SAaxis.posteriorcorrs) <- sprintf("draw%s",seq(from = 1, to = ncol(deriv.SAaxis.posteriorcorrs)))
deriv.SAaxis.posteriorcorrs <- cbind(deriv.SAaxis.posteriorcorrs, (gam.derivatives.glasser %>% group_by(age) %>% group_keys()))

Age of maximal sensorimotor-association axis correlation: posterior median value + 95% credible interval

age.max.corr <- deriv.SAaxis.posteriorcorrs %>% #find the age at which the correlation is largest for each draw
    summarise(across(contains("draw"),
                     .fns = function(x){
                       round(deriv.SAaxis.posteriorcorrs$age[which.max(x)],4)
                     }))
age.max.corr <- t(age.max.corr)
age.max.corr.median <- median(age.max.corr) #median age #bayes
age.max.corr.CI <- quantile(age.max.corr, probs = c(0.025, 0.975)) #compute the credible interval based on all draws
age.max.corr.lower <- age.max.corr.CI[[1]]
age.max.corr.upper <- age.max.corr.CI[[2]]
age.max.corr.median
## [1] 15.0243
age.max.corr.CI
##    2.5%   97.5% 
## 14.6516 15.3224
age.max.corr <- as.data.frame(age.max.corr)
ggplot(age.max.corr, aes(x = V1)) + 
geom_histogram(bins = 12, fill = "#FFE5B1") + 
geom_vline(xintercept = 15.0243, color = "#ba275f", size = .5) +
labs(x="\nAge", y="Posterior Draw Number\n") +
theme_classic() +
theme(axis.text.x = element_text(size = 6, family = "Arial", color = c("black")), axis.text.y = element_blank(), axis.title = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
scale_x_continuous(breaks=c(14.5, 15, 15.5))   

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/AgeMaxCorr_histogram.pdf", device = "pdf", dpi = 500, width = 1.03 , height = 1.02)

Age of first null (zero) sensorimotor-association axis correlation: posterior median value + 95% credible interval

age.zero.corr <- deriv.SAaxis.posteriorcorrs %>% #find the age at which the correlation first hits 0 for each draw
    summarise(across(contains("draw"),
                     .fns = function(x){
                      round(min(deriv.SAaxis.posteriorcorrs$age[round(x, digits = 1) == 0]),4)
                     }))
age.zero.corr <- t(age.zero.corr)
age.zero.corr.median <- median(age.zero.corr) #median age #bayes
age.zero.corr.CI <- quantile(age.zero.corr, probs = c(0.025, 0.975)) #compute the credible interval based on all draws
age.zero.corr.lower <- age.zero.corr.CI[[1]]
age.zero.corr.upper <- age.zero.corr.CI[[2]]
age.zero.corr.median
## [1] 19.273
age.zero.corr.CI
##    2.5%   97.5% 
## 18.6767 20.1675

Sensorimotor-association axis correlation value at each age: posterior median values + 95% credible interval

deriv.SAaxis.posteriorcorrs <- deriv.SAaxis.posteriorcorrs %>% select(contains("draw"))

deriv.SAaxis.mediancorr <- apply(deriv.SAaxis.posteriorcorrs, 1, function(x){median(x)}) #median correlation value at each age

cor.credible.intervals <- apply(deriv.SAaxis.posteriorcorrs, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws
cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals <- cbind(cor.credible.intervals, (gam.derivatives.glasser %>% group_by(age) %>% group_keys()))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr 
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")

Maximal correlation value

cor.credible.intervals %>% arrange(desc(SA.correlation)) %>% slice_head()
##       lower     upper      age SA.correlation
## 1 0.6629209 0.7040393 15.02429      0.6837556

Sensorimotor-association axis development sliding window plot

cor.credible.intervals$max.corr.CI <- (cor.credible.intervals$age > age.max.corr.lower & cor.credible.intervals$age < age.max.corr.upper) #add a column that indicates whether each age is included in the age of maximal correlation credible interval (T/F)
cor.credible.intervals$max.cor.window <- cor.credible.intervals$age*cor.credible.intervals$max.corr.CI #add a column that only includes ages in this interval
cor.credible.intervals$max.cor.window[cor.credible.intervals$max.cor.window == 0] <- NA
                        
cor.credible.intervals$zero.corr.CI <- (cor.credible.intervals$lower < 0 & cor.credible.intervals$upper > 0) #add a column that indicates whether the credible interval for the correlation includes 0
cor.credible.intervals$zero.corr.window <- cor.credible.intervals$age*cor.credible.intervals$zero.corr.CI #add a column that only includes ages in the zero window
cor.credible.intervals$zero.corr.window[cor.credible.intervals$zero.corr.window == 0] <- NA
ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
geom_line(size = .5) +
geom_ribbon(alpha = .2, fill = c("grey30")) +
labs(x="\nAge", y="Developmental Alignment to the Axis\n") +
geom_ribbon(aes(x = max.cor.window, y = SA.correlation), fill = "#ba275f") +
theme_classic() + 
theme(axis.text = element_text(size = 7, family = "Arial", color = c("black")), axis.title = element_text(size = 7, family = "Arial", color = c("black")), axis.line = element_line(size =.22), axis.ticks = element_line(size = .22)) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .1)) +
scale_y_continuous(breaks=c(0,0.2,0.4,0.6)) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure4/Derivative_SAaxis_correlation_plot.pdf", device = "pdf", dpi = 500, width = 3.82 , height = 2.35)

Individual differences in intrinsic fMRI fluctuations are associated with differences in youth’s developmental environment

Significant associations between fluctuation amplitude and environmental factors

Number of cortical regions showing significant environment effects

sum(p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05) 
## [1] 141

Percent of cortical regions showing significant environment effects

(sum(p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05))/(nrow(df.env))*100 
## [1] 41.96429

Across-cortex environment effect size and direction

ggseg(.data = df.env, atlas = "glasser", mapping=aes(fill=GAM.env.tvalue), position = c("stacked")) + 
  theme_void() + 
  labs(fill="Age Effect") +
  scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits=c(-5,5), oob=squish) +
  theme(legend.text = element_text(size = 6, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure6/Brainmap_enveffect.jpg", device = "pdf", dpi = 500, width = 3.05 , height = 2.25)

Correlation between environmental effects and the sensorimotor-association axis

#all regions
cor.test(df.env$GAM.env.tvalue, df.env$SA.axis, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  df.env$GAM.env.tvalue and df.env$SA.axis
## S = 3274572, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4820453
perm.sphere.p(df.env.spin$GAM.env.tvalue, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')
## [1] 0
#significant regions
df.env.spin$GAM.env.tvalue.significant <- df.env.spin$GAM.env.tvalue*(p.adjust(df.env.spin$Anova.env.pvalue, method=c("fdr")) < 0.05)
df.env.spin$GAM.env.tvalue.significant[df.env.spin$GAM.env.tvalue.significant == 0] <- NA

cor.test(df.env.spin$GAM.env.tvalue.significant, df.env.spin$SA.axis, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  df.env.spin$GAM.env.tvalue.significant and df.env.spin$SA.axis
## S = 204474, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5623229
perm.sphere.p(df.env.spin$GAM.env.tvalue.significant, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')
## [1] 0

Correlation plot

df.env$significant <- p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05

ggplot(df.env, aes(x = SA.axis, y = GAM.env.tvalue, fill = GAM.env.tvalue, color = significant)) + 
geom_point(size = 1.5, stroke = .12, pch = 21) +
scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits = c(-5,5), oob = squish) +
scale_color_manual(values = c("white","black")) +
labs(x="\nSensorimtor-Association Axis Rank", y="Environment Effect\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .5) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22)) +
scale_y_continuous(breaks = c(-5, 0, 5))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/Neuron/Figure6/AxisEnv_correlationplot.jpg", device = "pdf", dpi = 500, width = 2.5 , height = 1.7)